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Abstract. In this paper we discuss recent applications of the Smoothed Par- 
ticle Hydrodynamics (SPH) method to the simulation of supersonic turbulence 
in the interstellar medium, as well as giving an update on recent algorithmic de- 
velopments in solving the equations of magnetohydrodynamics (MHD) in SPH. 
Using high resolution calculations (up to 134 million particles), we find excellent 
agreement with grid-based results on a range of measures including the power 
spectrum slope in both the velocity field and the density-weighted velocity p 1 ^ 3 v, 
the latter showing a Kolmogorov-like fc~ 5 / 3 scaling as proposed by Kritsuk et 
al. (2007). We also find good agreement on the statistics of the Probability 
Distribution Functio n (PDF) and structure function s , inde pendently confirming 
the scaling found bv lSchmidt. Federrath fc'K lesscn (2008). On Smoothed Par- 
ticle Magnetohydrodynamics (SPMHD) we have recently wasted a great deal of 
time and effort investigating the vector potential as an alternative to the Euler 
potentials formulation, in the end concluding that using the vector potential has 
even more severe problems than the standard (B-ficld based) SPMHD approach. 



1. Introduction 

Turbulence and magnetic fields are thought to be two of the most important 
ingredients in the star formation process, so much so that there remains ongoing 
debate — both observational and theoretical — as to which one is the controlling 
factor (e.g. lMac Low & Klessenll2004l ; iKunz k Mouschovias!l2009l ) . We have only 



recently begun to understand the role of either in depth, primarily as a result 
of our increased ability to simulate both in numerical calculations of the star 
formation process. 

Several authors have proposed that turbulence in the interstellar medium 
(ISM) can be used to predict the form of the st ellar initial mass function 
(IMF) thus providing a 'the ory' of star formation (jPadoan fc Nordlundl [2002 ; 



iHennebelle Sz Chabrierl 12008). Any such theory requires consensus on the ba- 
sic statistical characteristics of turbulence in the supersonic regime and as- 
sumes that these are universal — that is, independent of boundary condi- 
tions and_djivmg^jneclianisms — for which there is some observational support 
(e.g. Hever Sz Brunt 20041 ). However there exists considerable disagreement be- 
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tween results obtain ed using different numerical codes, most recently between 
iPadoan et al. ( 2007) who claim tha t the statistics of turbulence are universal 
and iBallesteros-Paredes et al. ( 20061 ) who claim that they are not, based on cal- 
culations utilising both a Smoothed Particle Hydrodynamics (SPH) code and a 
grid-based Total- Variation-Diminishing (TVD) code. 

This kind of disagreement, and the need within star formation theory for a 
consensus on the basic statistics of supersonic turbulence has prompted at least 
two maj or code comparison pr ojects over the last few years, the "Potsdam" com- 
parison (jKitsionas et al.1 12009). comparing decaying, hydrodynamic turbulence 
and the KITP comparison (unpublished), comparing decaying hydrodynamic 
and magnetohydrodynamic (MHD) turbulence. However these comparisons suf- 
fer from the limited time evolution that can be obtained from a decaying tur- 
bulence simulation as well as the numerical issues posed by starting from an 
evolved snapshot produced by a particular code. 

We have therefore undertaken our own very detailed comparison of driven 
turbulence using just two codes, an SPH code, phantom, and a grid code, FLASH 
(used in uniform grid mode) taken to be broadly representative of their class of 
codes, in order to see whether agreement can be reached (and at what resolution) 
on the statistics of su personic turbulence appr opriate to the ISM. The results 
are published in full in lPrice &: Federrath ( 2010l ) and we only provide a snapshot 
here of the main results (^EQ), referring the reader to that paper for more detailed 
information. 

What we are not yet able to address is the combination of driven turbu- 
lence and MHD in SPH, primarily because of the limitations posed by the Euler 
potentials formulation used as the only method that sufficiently maintains the 
divergence-free (V-B = 0) constraint o n the magnetic field such that st ar forma- 
tion calculations can be performed (e.g. iPrice fc Batdl2007l . 120081 . [2009) without 
the stars themselves exploding due to numerical errors (see ffiT]). Following a 
suggestion from Axel Brandenburg, we therefore embarked on a quest to ex- 
amine whether or not the use of the vector potential could similarly solve the 
divergence problem in the context of Smoothed Particle Magnetohydrodynamics 
(SP MHD) without the associated topological restrictions of the Euler potentials 
(see iBrandenburd |2010D . In short, this turned out not to be the case (read the 
horror that is |Pricdl2010l ). However, for the pleasure of the reader we describe 
the tortuous journey in 



2. Turbulence: Grids vs. SPH 

We adopted the usual approach to the artificial production of turbulence in 
numerical codes, using periodic boundary conditions and applying a "large- 
scale" driving in fourier space with random amplitudes and phases that are 
slowl y changed over a tim escale roughly comparable to the box crossing time 
(see iFederrath et al.l 12008 for details). In order to obtain as close a comparison 



as possible we pre-generated the random driving pattern for all times during 
the simulation such that both codes use an identical forcing. Studying driven 
turbulence means that there are no numerical issues in either code relating to 
the (very simple) uniform initial conditions. 
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Figure 1. Column density in Mach 10 turbulence calculations employing 
512 3 computational elements using the SPH code (phantom, left), using the 
grid code (flash, middle), and using an SPH density calculated from the 
FLASH tracer particles (right). 



Figure [T] shows the column density after two dynamical times (adopting the 
usual definition = L/(2A4) where L is the box size and M is the RMS Mach 
number) from Mach 10 calculations employing 512 3 computational elements in 
both grid cells and SPH particles. Whilst 512 3 is only a moderate resolution in 
current grid-based turbulence simulations, our use of 512 3 or approximately 134 
million SPH particles makes this easily the highest resolution turbulence calcula- 
tion ever performed with SPH. Indeed, given that SPH is a much more expensive 
comp utational method compared to uniform-grid calculations (iKitsionas et al.l 
20091 ) it is reassuring that dense regions appear significantly better resolved in 



the SPH results (with adaptive smoothing lengths the effective SPH resolution 
at any point is ~ re 1 ' 3 /!/ where re is the local particle density). Since any 
viscosity in SPH must be explicitly added to the calculation (in the form of ar- 
tificial viscosity, AV), the Reynolds number at any point in the flow can also be 
explicitly calculated. Figure [2] shows the line-of-sight averaged Reynolds num- 
ber, (Re) = /(pRe)dz/ J p dz, where the Reynolds number for each particle a 
has been calculated according to Re = \v\L/is, where \v\ is the magnitude of the 
particle vel ocity, L is the box size and the viscos ity v is given by v = 1/WotAvCsh 
in 3D (e.g. lMurravlll996l ; lLodato fe PricdfeoiOl ). The mean Re is around 10 5 in 
the highest resolution calculation, with individual Re's reaching ~ 10 6 in the 
densest regions. 

Another unique aspect of our comparison was adding Lagrangian tracer 
particles to the flash calculation, where we have used the SPH density calcu- 
lation routine from phantom in order to calculate a density field for the tracer 
particles distribution (right panel of Figure [T]) . Doing so reveals an amazing 
resolution of sub-grid structures compared to the raw grid-based density field 
(middle panel). The caveat is that in the Probability Distribution Function it is 
clear that at least some of this structure may be an artefact of the tracer particle 
dynamics below the grid scale, in that the density PDF on the tracers in the low 
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Figure 2. Line-of-sight averaged Reynolds number J(pRe)dz/ J pdz in the 
512 3 SPH Mach 10 turbulence calculation. The mean Reynolds number at this 
resolution is around 10 5 , with individual Reynolds numbers reaching ~ 10 6 
in the densest regions. 
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Figure 3. Compensated, time averaged power spectra in the 512 3 grid 
(dashed, red) and SPH (solid, black) turbulence calculations, showing the 
velocity field (top, compensated by k 2 ) and for the quantity p 1 / 3 ?; (bottom, 
compensated by /c 5 / 3 ) which appears to scale in a Kohnogorov-like way. 
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resolution runs does not appear to match the grid-based density PDF at higher 
resolution. 

Whilst the SPH is generally better reso lved in the density field, reflected 
also in the PDFs (see Price &: Federrath 20101 ) , it is not necessarily true for more 
volumetric statistics, such as power spectra (Figure [3]) where we find comparable 
results roughly when the number of computational elements is equal (Figure EJ. 
The power spectra of the density-weighted quantity p l ^v show a Kolmogorov- 
like /c -5 / 3 s caling at high resolu tion, in agreement with the 'new universality' 
proposed bv lKritsuk et all (|2007l 1. SPH appears more dissipative in pure velocity 
statistics (top panel of Figure [3]) but correspondingly less dissipative in density- 
weighted statistics (bottom panel of Figure [3]) . 

Overall, we find very good agreement between the SPH and grid results 
on the statistics of supersonic turbulence, provided sufficient resolution is em- 
ployed in either code, where "sufficient" depends on the desired statistic — we 
require roughly 512 3 computational elements in either code in order to resolve 
the inertial range in the power spectrum, whereas a good estimate of the PDF 
is obtained at lower resolutions. We find excellent agreement at high resolution 
be tween our measured structure functio n scaling and results obtained recently 
bv lSchmidt. Federrath fc Klessenl (120081 ) at Mach 6. 



3. Smoothed Particle Magnet ohydro dynamics: an update 

Whilst SPH is widely used for solving the equations of hydrodynamics in many 
astrophysical applications, the incorporation of magnetic fields into the SPH 
algorithm has a long and somewhat tortuous history. This is primarily due to a 
numerical instability that occurs if one tries to formulate the SPMHD equations 
such that momentum is conserved exactly (i.e., as the gradient of a stress tensor). 
In particular, SPH in its most basic form relies on a net repulsive force between 
particles to preserve order in the particle distribution. In MHD the stress tensor 
can become negative when the magnetic pressure is greater than the gas pressure 
so the particles can attract each other unstoppably. However physically any 
such force along the line of sight is spurious, related to the non-zero divergence 
of the magnetic field. Amongst other things (e.g. formulating shock-capturing 
dissipative terms, incorporating spatially variable smoothing lengths), dealing 
with this "cl umping instability" was one o f the key steps to a robust algorithm 
developed in lPrice fe Monaghanl (l2004al lbl. 120051 ). 



A related but not identical issue, particularly for 3D calculations, is enforce- 
ment of the divergence- free constraint. Already what one means by "divergence- 
free" in SPH is unclear, because the stencil for measuring divergence of a vector 
field depends not only on the operator chosen but also on the particle distri- 
bution which can be different from one timestep to the next. It is primarily 
for this reason that the d ivergence cleaning schemes such as those discussed in 
Price Sz Monaghanl (2005) are not very successful for 3D problems. Far bet- 



ter are approaches that have the divergence constraint "built-in". In particular 
using the 'Euler potentials' formulation, B = Va x V/3, is a very elegant ap- 
proach for a Lagrangian code, not only because the divergence is zero by con- 
struction but also because the induction equation adopts the very simple form 
da/dt = 0;d/3/dt = corresponding to the advection of field lines by Lagrangian 
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Figure 4. The meaning of "maintaining the divergence constraint" in 
Smoothed Particle Magnetohydrodynamics. In the absence of effective control 
of divergence errors (left) the collapsing protostar explodes due to large nu- 
merical errors in the force in the region of extreme magnetic field gradients. 
Using the Euler potentials (right) these errors can be effectively controlled 
(and simulations followed beyond the point of star formation), though physi- 
cal winding of the field is also lost. 



particles. Use of the Euler potentials has enabled us to study the role o f mag- 
netic fields in realistic 3D star formation problems ([Price Bate! 120071 . 120081 . 



120091 ). where the meaning of not maintaining the divergence constraint is clear. 
Figure H] gives an example of this, showing the results of a spherical, rotating 
collapse calculation with an initially uniform ma gnetic field threading the i nitial 
cloud. Using a "standard" SPMHD approach ([Price &: Monaghanl 120051 ). the 
protostar simply explodes due to numerical force errors caused by the extremely 
large gradients and associated non-zero divergence in the magnetic field (left 
panel). Using the Euler potentials one is able to follow the calculation stably to 
well beyond the point of star formation (right panel). 

However, the "advection" property of the Euler potentials also means that 
physical winding of the fields is lost once a one-to-one mapping of the field from 
the initial conditions ceases to exist. For this reason we have recently explored an 
alternative approach based on the vector potential, the hope being that similar 
maintenance of the divergence-free condition could be achieved without the loss 
of physical winding processes. However the induction equation for the vector 
potential takes a rather more complicated form: 

dA 

— = v xVxA + (v- V) A -T]J + V0, (1) 

(JUL 

where v is the velocity, r] is the resistivity, J is the current and <f> is an arbi- 
trary scalar corresponding to the freedom to choose a gauge. The second term is 
particularly nasty since it represents effectively trying to compute "reverse ad- 
vection" on the particles — the absence of an advective term being the primary 
advantage of Lagrangian schemes over their Eulerian counterparts. Nevertheless 
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Axel Brandenburg suggested that choosing eft = A • v as a gauge might solve 
these problems since the gradients are flipped onto the v instead of A. Choosing 
such a gauge also means that the induction equation becomes Galilean invari- 
ant, which in turn implies that it can be used to derive an exactly-conservative 
SPMHD formulation based on the vector potential ab initio using a variational 
principle. 

The great hope was that such a vector potential-based formulation, with the 
divergence-constraint "built-in" as in the Euler potentials method, could provide 
an elegant method which was both conservative and stable, without any of the 
numerical problems of previous formulations. Alas, this turned out not to be the 
case — in fact the equations so derived (see |PriceIl201Cll ). whilst elegant, turned 
out to be even more unstable to the clumping instability than in the standard 
approach. Furthermore it was found that another, independent instability occurs 
related to the unconstrained growth of unphysical vector potential components, 
most likely related to poor accuracy with respect to the "reverse- ad vection"- 
type terms. The latter means that the vector potential cannot be used even in 
conjunction with a standard (stable but non-conservative) SPMHD force term. 
The somewhat painful conclusion was therefore that the vector potential is not 
a viable approach for SPMHD. 
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